


% Iterate on lower boundary:
dJ_tol   = 1e-9;
dJ_error = 1;
biter    = 1;
maxbiters= 50;
while dJ_error > dJ_tol && biter <=maxbiters
    
    mL = .5*(mLlo + mLhi);      % [mLlo, mLhi] from solve_initialize.m
    

    %% Natural wastage
     % Solve for J1, J2 from value matching (VM) and smooth pasting (SP) conditions @ mL, given m_L
    muhat  = mu + (1-alpha)*s*lambda;
    rho0   = r  + s*lambda;
    rho1   = rho0 - muhat;
    gamma1 = ( -(muhat-((sigma^2)/2)) - sqrt( (muhat-((sigma^2)/2))^2 + 2*(sigma^2)*rho0 ) )/(sigma^2);
    gamma2 = ( -(muhat-((sigma^2)/2)) + sqrt( (muhat-((sigma^2)/2))^2 + 2*(sigma^2)*rho0 ) )/(sigma^2);
    J2     = ((omega0/rho0) - ((1-omega1)/rho1)*((gamma1-1)/gamma1)*mL) * (gamma1/(gamma1-gamma2)) * (mL^-gamma2);
    J1     = -((mL^(1-gamma1))/gamma1) * ( ((1-omega1)/rho1) + gamma2*J2*(mL^(gamma2-1)) );
     % Next, solve VM condition @ mM boundary
    x  = linspace(mL+.01, mL*2.5, 100)';
    residm = c - ((1-omega1)/rho1)*x + (omega0/rho0) - J1*(x.^gamma1) - J2*(x.^gamma2);
    if min(residm) >0
        disp('ERROR: No solution to value matching problem @ mM (boundaries)');
        pause;
    else
        [~,pos] = min(residm);
        if pos >1 || pos <length(x)
            % interior min
            mhat = x(pos);
        elseif pos==length(x)
            mhat = x(end);
        elseif pos==1
            disp('ERROR: No solution to value matching problem @ mM (boundaries)');
        end
    end
     %
     % mM must between mL and mhat --> Define x s.t. mM = mL + (mhat-mL)*(exp(x)/(1+exp(x)))
     % Initial guess: mM = 1.25*mL. Use this to form initial guess for x and solve for x.
    mM0 = 1.25*mL;
    lhs = (mM0 - mL)/(mhat-mL);
    x0  = log(lhs/(1-lhs));
    [xM, fMval, fMflag] = fzero(@(x) (c - ((1-omega1)/rho1)*(mL+(mhat-mL)*(exp(x)/(1+exp(x)))) + (omega0/rho0) - J1*((mL+(mhat-mL)*(exp(x)/(1+exp(x)))).^gamma1) - J2*((mL+(mhat-mL)*(exp(x)/(1+exp(x)))).^gamma2)), x0);
    mM = mL + (mhat-mL)*(exp(xM)/(1+exp(xM)));
    if fMflag <1
        disp('ERROR: Failed to solve for mM boundary (boundaries)');
        pause;
    end 
    g  = mM/mL;
     % For later, Evaluate J'(m_M) in natural wastage region (left side derivative).
    dJ_mM_left  = ((1-omega1)/rho1) + gamma1*J1*(mM^(gamma1-1)) + gamma2*J2*(mM^(gamma2-1));



    %% Full Replacement
     % First, solve JJ1, JJ2 using VM equs for mM and mR, J(mM) = c and J(mR) = c+C
     % mM is already given. Solve over grid of mR. Then, solve SP equ. for mR
    Rho0 = r;
    Rho1 = r-mu;
    phi1 = ( (.5*(sigma^2)-mu) - sqrt( (.5*(sigma^2)-mu)^2 + 2*(sigma^2)*r) ) /(sigma^2);
    phi2 = ( (.5*(sigma^2)-mu) + sqrt( (.5*(sigma^2)-mu)^2 + 2*(sigma^2)*r) ) /(sigma^2);
    scriptw = ((1/(1-alpha)) - phi1) / (phi2-phi1);
     %
    aux1  = ((1-alpha)/((sigma^2)/2));
    mRobs = 5000;	% Add more grid points at this stage to increase accuracy of mR solution.
    quadobs = 50;   
    mRvec = linspace( mM+1e-4, mRmax, mRobs)';      % mRmax is from solve_initialize.m
    mRmat = ones(quadobs,1)*mRvec';
    Gvec  = mRvec/mM;
    Gmat  = mRmat/mM; 
     %
    [nodes,wghts] = legendre(mM, mRvec, quadobs);   % quadobs by mRobs
    delta    = s*lambda ./ (1 + aux1*s*lambda*log(nodes/mM));
    integrand=  (  scriptw *((mRmat./nodes).^phi1) + ...
                (1-scriptw)*((mRmat./nodes).^phi2) ).*(delta./nodes);
    R_mR = c*aux1 * diag(wghts'*integrand);         % replacement costs
    JJ1  = C +(c+(omega0/Rho0))*(1-Gvec.^phi2) - ((1-omega1)/Rho1)*(Gvec - Gvec.^phi2)*mM + R_mR;
    JJ1  = JJ1 ./ ( (mM^phi1)*(Gvec.^phi1 - Gvec.^phi2) );
    JJ2  = (c - ((1-omega1)/Rho1)*mM + (omega0/Rho0) - JJ1*(mM^phi1)) * (mM^-phi2);
     %
     % Now solve for mR using SP condition
    integrand = (phi1*   scriptw *((mRmat.^(phi1-1))./(nodes.^phi1)) + ...
                 phi2*(1-scriptw)*((mRmat.^(phi2-1))./(nodes.^phi2))).*(delta./nodes);
    dR_mR     = c*aux1*(  (delta(end,:)'./mRvec) + diag(wghts'*integrand)  ); 
    resid  = ((1-omega1)/Rho1) - dR_mR + phi1*JJ1.*(mRvec.^(phi1-1)) + phi2*JJ2.*(mRvec.^(phi2-1)); 
    if resid(end) >0
        disp('ERROR: Unable to solve for mR boundary (boundaries)');
        pause;
    else
        mR     = interp1(resid,mRvec,0);
        G      = mR/mM;
        JJ1    = interp1(mRvec, JJ1, mR);
        JJ2    = interp1(mRvec, JJ2, mR);
    end
     %
     % Evaluate dJ/dm @ mM in full replacement region (right side derivative):
    dR_mM = c*aux1 *s*lambda/mM;
    dJ_mM_right = ((1-omega1)/Rho1) - dR_mM  + phi1*JJ1*(mM^(phi1-1)) + phi2*JJ2*(mM^(phi2-1));
    dJ_diff     = dJ_mM_left - dJ_mM_right;  
     %
     % Update guess for mL until left- and right-side derivatives of J @ mM equate. 
    if dJ_diff <0 
        % dJ_diff is decreasing in mL -->
        % mL too large
        mLhi = mL;
    else
        mLlo = mL;
    end
    dJ_error = abs(dJ_diff);
    biter = biter +1;
end
 %
 %
if biter < maxbiters
     % Upper boundary in net expansion region, delta(mU) =0
    delta_mR   = s*lambda ./ (1 + aux1*s*lambda*log(mR/mM));
    const      = ((1-omega1)*mR/alpha) - omega0 - r*C - (r+delta_mR)*c;
    [mU,fUval, exitflag] = fzero(@(x) (delta_mR + (1/c)*( (((1-omega1)*(x - mR))/alpha) - const*( (x/mR).^(1/(1-alpha)) - 1 ) )), mR*1.1 ); 
    if exitflag <1
        disp('ERROR: Failed to solve for mU');
        pause;
    end
    if mU <mR
        disp('ERROR: mU < mR');
    end
     
else
    disp('ERROR: Unable to solve for boundaries');
    disp('Residual and tolerance='); 
    disp( [dJ_error,dJ_tol] ); 
    pause;
end


